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ABSTRACT 

The merger of binary neutron stars is likely to lead to differentially rotating remnants. In this paper we 
numerically construct models of differentially rotating neutron stars in general relativity and determine 
their maximum allowed mass. We model the stars adopting a polytropic equation of state and tabulate 
maximum allowed masses as a function of differential rotation and stiffness of the equation of state. We 
also provide a crude argument that yields a qualitative estimate of the effect of stiffness and differential 
rotation on the maximum allowed mass. 
Subject headings: Gravitation — relativity — stars: rotation 



1. INTRODUCTION 

Neutron stars that are newly formed in the coalescence 
of binary neutron stars (or in supernova collapse) are likely 
to be differentially rotating (see, e.g., the fully dynamical 
simulations of Rasio & Shapiro 1992,1994, Shibata & Uryu 
2000, Faber, Rasio & Manor 2001, as well as the review of 
Rasio & Shapiro 1999 and references therein). Differential 
rotation may play an important role for the stability of 
these remnants, since it can be very effective in increasing 
their maximum allowed mass. 

Assuming that neutron stars in binaries have individual 
masses close to 1.4 Mq (Thorsett & Chakrabarty 1999), 
and assuming that the maximum allowed mass of a non- 
rotating neutron star is in the range of 1.8 - 2.3 Mq (Ak- 
mal, Pandharipande & Ravenhall 1998), one might con- 
clude that the coalescence of binary neutron stars leads to 
immediate collapse to a black hole. However, both thermal 
pressure and rotation can increase the maximum allowed 
mass. 

For coalescence from the innermost stable circular or- 
bit, thermal pressure is believed to have a small effect. 
The maximum mass of uniformly rotating stars is lim- 
ited by the spin rate at which the fluid at the equator 
moves on a geodesic; any further speed-up would lead to 
mass shedding (the "Kepler" limit). Uniform rotation can 
therefore increase the maximum allowed mass by at most 
about 20 % for very stiff equations of state (Cook, Shapiro 
& Teukolsky, 1992, hereafter CSTl, and 1994, hereafter 
CST2), which is not sufficient to stabilize remnants of bi- 
nary neutron star merger. Uniformly rotating equilibrium 
configurations with rest masses exceeding the maximum 
rest mass of nonrotating stars constructed with the same 
equation of state are referred to as "supramassive" stars 
(CSTl, CST2). 

Differential rotation, however, can be much more effi- 
cient in increasing the maximum allowed mass. In differ- 
entially rotating stars, the core may rotate faster than the 
envelope, so that the core can be supported by rapid ro- 
tation without the equator having to exceed the Kepler 
limit. This effect was demonstrated in Newtonian grav- 
itation by Ostriker, Bodenheimer & Lynden-Bell (1966) 



for white dwarfs, and in general relativity by Baumgarte, 
Shapiro & Shibata (2000, hereafter BSS) forn = 1 poly- 
tropes. BSS also showed by way of example that stars with 
about 60 % more mass than the maximum allowed mass 
of the corresponding non-rotating star can be dynamically 
stable against both radial and nonaxisymmetric modes. 
BSS refer to differentially rotating equilibrium configura- 
tions with rest masses exceeding the maximum rest mass 
of a uniformly rotating star as "hypermassive" stars. 

In this paper we extend the findings of BSS for n — I 
polytropes in two ways. We survey different polytropic in- 
dices, and study the effect of both differential rotation and 
stiffness of the equation of state on the maximum allowed 
mass. We also introduce a very simple model calculation 
that illustrates these effects qualitatively, and that more- 
over gives surprisingly accurate results for soft equations 
of state and moderate degrees of differential rotation. 

This paper is organized as follows. In Section 2 we 
present qualitative considerations, leading to a simple esti- 
mate of the maximum mass of differentially rotating neu- 
tron stars. In Section 3 we construct numerical models of 
fully relativistic, differentially rotating neutron stars for 
different polytropic indices. We briefly discuss our find- 
ings in Section 4. We also include Appendix A with tables 
of our numerical results. 

2. QUALITATIVE CONSIDERATIONS 

When magnetic fields and relativistic effects can be ne- 
glected, rotating stars in equilibrium satisfy the Newtonian 
virial theorem 

PF + 2T + 3n = Vu(^l-2|^j +3H = (1) 

(see, e.g., Shapiro & Teukolsky 1983). Here the rotational 
kinetic energy T scales as 

T~^, (2) 



MRr 



the potential energy as 



W 



Rr 



(3) 



2 



LYFORD, BAUMGARTE & SHAPIRO 



and the internal energy, computed from a volume integral 
of the pressure, as 

n ~ p^/"M. (4) 

where J is the angular momentum, M the mass, the 
equatorial radius and p the mass density. We have also 
assumed a polytropic equation of state 

P = Kp^+^'''\ (5) 

where K ]s a. constant and n the polytropic index. Here 
and throughout we set the gravitation constant G =1. 

Inserting the above expressions into the virial equation 
(1), we find 



(6) 



where [3 and 7 arc the appropriate coefficients. In general, 
/3 depends on the eccentricity of the star, but restricting 
our analysis to small values of T/|VF|, and hence to nearly 

spherical stars, we assume that both (3 and 7 arc constant. 

For non-rotating stars (T = 0) we solve equation (6) to 
find 

3/2 



M = 



(3-n)/(2n) 



while for rotating stars we obtain 



M.nt, = M 1 - 2 



T 

'W\ 



-3/2 



(7) 



(8) 



For small values of T/|M^| the right hand side can be ex- 
panded, and we then find for the mass increase 6M = 
Mrot - M 

(cf. Shapiro & Teukolsky 1983, eqn 7.4.40). This expres- 
sion determines the fractional mass increase as a function 
of T/lVFj for a constant value of the density p. We will 
use this result to estimate the increase in the maximum 
allowed mass of a neutron star, even though typically ro- 
tating stars assume their maximum masses at slightly dif- 
ferent densities than the corresponding non-rotating stars. 

To evaluate T and \W\ we further simplify the prob- 
lem by assuming that the star's density profile is a step 
function with a constant density pc (equal to the original 
central density) inside a spherical core of radius i?c, and 
zero outside. Prom requiring that the mass of this model 
star, 

47r 

-PcRl (10) 



M = 



be equal to the original mass 
M = — 



(11) 



where p is the average density, we find the following rela- 
tion between the central condensation and the ratio of the 
radii 

Re ^ (Pc^ 
Rc \P 

Further assuming that the core is uniformly rotating with 
the central angular velocity fic, we find 



(12) 



-.m 



MRt 



(13) 



Inserting this relation together with the potential energy 

3M2 



W 



5 Rr. 



into (9) yields 



5M Rl^l 
M ^ M ' 



(14) 



(15) 



To find the increase in the maximum allowed mass, it is 
useful to assume that the star rotates at the mass-shedding 
limit 



M 

Rr 



(16) 



where fig is the equatorial angular velocity. This equation 
can now be used to eliminate the mass M in equation (15), 
which, together with (12), yields 



6M 



p_ fn, 

Pc 



(17) 



Equation (17) provides a very simple estimate for the in- 
crease of the maximum allowed mass. It depends only on 
the central condensation of the non-rotating star, which is 
a function of the stiffness of the equation of state, and the 
ratio of the angular velocities at the center and equator, 
which is a function of the the degree of differential rota- 
tion. For uniformly rotating stars, the maximum mass 
increase is estimated to be simply the inverse of the cen- 
tral condensation. In Table 1 we compare this estimate 
with the numerical findings of CST2 and find remarkably 
good agreement for soft equations of state. Table 1 also 
illustrates an ambiguity; in Newtonian gravity, the central 
condensation is uniquely determined by the polytropic in- 
dex, but in general relativity the central condensation of 
a star depends on the central density. We therefore com- 
pute a "relativistic" central condensation Pc/p|tov from 
the central energy density Pc and an average density de- 
fined as 



P = 



3M 
47ri?3 



(18) 



of the non-rotating maximum mass model, where M is the 
total mass-energy of the star, and R is the circumferential 

radius. We find that this value yields better agreement 
with the numerical values of the maximum mass increase 
than adopting the Newtonian central condensation. 

The ratio T/|iy | provides a useful criterion for the onset 
of secular (T/|W| ~ 0.14) or dynamical (T/|VF| ~ 0.27) 
non-axisymmetric instabilities. Inserting these limits into 
equation (8) shows that mass increases of secularly stable 
stars are limited by 5M/M < 1.63 while mass increases 
of dynamically stable stars are limited by SM/M < 3.2. 
These values agree quite well with the respective limits of 
1.70 and 3.51 found by Shapiro & Teukolsky (1983; equa- 
tions (7.4.41) and (7.4.42)), who also take stellar deforma- 
tions into account. 
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Table 1 

Maximum Mass increase for uniformly rotating polytropes. 





^max b 


j^max c 


d 

f'max 


SM/M\csT " 


p/Pc\tov * 


P/Pc Newt ^ 


0.5 


0.151 


0.395 


1.29 


0.224 


0.375 


0.545 


1.0 


0.180 


0.763 


0.42 


0.146 


0.209 


0.304 


1.5 


0.276 


1.97 


0.072 


0.099 


0.115 


0.167 


2.0 


0.523 


6.94 


5.8x10-3 


0.066 


0.063 


0.088 


2.5 


1.25 


41.0 


1.26x10-"* 


0.040 


0.034 


0.043 


2.9 


3.23 


620 


1.54x10-'' 


0.023 


0.021 


0.021 



"Polytropic index. 

''Maximum rest mass of non-rotating polytrope (CST2). 

^Circumferential radius of the non-rotating maximum-mass configuration (CST2). 

'^Maximum energy density of the non- rotating maximum mass configuration (CST2) . 

'^Fractional rest mass increase (CST2). 

■''Estimate (17) using relativistic central condensation. 

^Estimate (17) using Newtonian central condensation. 



3. NUMERICAL RESULTS 

We use a modified version of the numerical code of CSTl 
and CST2 to construct models of differentially rotating 
neutron stars. The code is based on similar algorithms 
developed by Hachisu (1986) and Komatsu, Eriguchi and 
Hachisu (1989), and we refer to CSTl for details. We 
adopt a polytropic equation of state 

P = Kp'+'/\ (19) 

where po is the rest-mass density and where equation (19) 
reduces to equation (5) in the Newtonian limit. We take 
the polytropic constant K to be unity without loss of gen- 
erality. Since iiT"/^ has units of length, all solutions scale 
according to M = if-"/^M, po = K"po, etc., where the 
barred quantities are dimensionless quantities correspond- 
ing to K = 1, and the unbarred quantities are physical 
quantities (compare CSTl). 

Constructing differentially rotating neutron star models 
requires choosing a rotating law F{^1) = u^Urf,, where u* 
and are components of the four-velocity u" and fl is 
the angular velocity. We follow CSTl and assume the ro- 
tation law F{U) = A'^{Uc — ri), where the parameter A 
has units of length. Expressing and in terms and 
metric potentials yields equation (42) in CSTl, or, in the 
Newtonian limit. 



l + i-2f2 sinks' 
Here we have rescaled A and r in terms of the equatorial 
radius Re- A = A/ Re and f = r/Re- The parameter A is 
a measure of the degree of differential rotation and deter- 
mines the length scale over which J7 changes. Since uni- 
form rotation is recovered in the limit ^ ^ oo, it is conve- 
nient to parametrize sequences by A~^ . In the Newtonian 
limit, the ratio ilc/^^e that appears in the estimate (17) is 
related to A~^ by Oc/^e = 1 + A~'^, but for relativistic 
configurations this relation holds only approximately. 

We adopt this particular rotation law for convenience 
and for easy comparison with many other authors who 
have assumed the same law. We also compared with 
the remnants' angular momentum distribution in the fully 



relativistic dynamical merger simulations of Shibata and 
Uryu (2000) and to the post-Newtonian simulations of 
Faber, Rasio and Manor (2001). We have found that their 
numerical results can be fit reasonably well by our adopted 
differential rotation law. 

We modify the numerical algorithm of CSTl by fixing 
the maximum interior density instead of the central den- 
sity for each model. This change allows us to construct 
higher mass models in some cases, since the central den- 
sity docs not always coincide with the maximum density 
and hence may not specify a model uniquely. 

For a given a value of n and A, we construct a se- 
quence of models for each value of the maximum density 
by starting with a static, spherically symmetric star and 
then decreasing the ratio of the polar to equatorial radius, 
Rpe = Rp/Re, in decrements of 0.025. This sequence ends 
when we reach mass shedding (for largo values of A), or 
when the code fails to converge (indic;atiiig the termina- 
tion of equilibrium solutions) or when Rpe = (beyond 
which the star would become a toroid). For each one of 
these sequences the maximum achieved mass is recorded. 
We repeat this procedure for different values of the max- 
imum density, covering about a decade below the central 
density of the non-rotating maximum mass model, which 
yields the maximum allowed mass for the chosen values of 
n and A. Our numerical results are tabulated in Appendix 
A. Our maximum mass increases are lower limits in the 
sense that even higher mass models may exist, but that 
we have not been able to construct them numerically. 

4. discussion and summary 

Our numerical results are tabulated in Tables A2 to A9 
in Appendix A. We also compare the increases in the 
maximum allowed mass with the estimate (17), and find 
surprisingly good agreement for soft equations of state and 
moderate degrees of differential rotation. 

In particular, we find that the fractional maximum rest 
mass increase 6M/M for uniformly rotating stars is well 
approximated by the inverse of the central concentration 
(see also Table 1). For moderate degrees of differential 
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rotation, 6M/M increases approximately with the square 
of the ratio between the central and equatorial angular 
velocity Qc/^e, in accord with equation (17). 

For all equations of state we find that 6M/M increases 
with ilc/^e only up to a moderate value of Qc/^e, and 
starts to decrease again for larger values (at least with our 
code and algorithm we do not find monotonically increas- 
ing mass configurations). For stiff equations of state this 
turn-around occurs for smaller values of llc/^e than for 
soft equations of state. For moderate degrees of differen- 
tial rotation ^c/^e ^ 2, a given value of fic/^^e will lead to 
a larger increase in the maximum allowed mass for a stiffer 
equation of state, as expected from the estimate (17). 

We find the largest maximum mass increases for moder- 
ately stiff equations of state. Some of these configurations 
exceed the maximum allowed mass of the corresponding 
non-rotating star by more than a factor of two. These con- 
figurations typically have large values of T/|VK| > 0.27, 
indicating that such stars may by dynamically unstable 
against bar formation (but see Shibata, Karino & Eriguchi 
2002, who found mild bar mode instabilities at very small 
values of for extreme degress of diff'erential rota- 

tion). They are also "toriodal", i.e. assume their max- 
imum density on a torus around the center of the star, 
which may indicate an m = 1 instability at even smaller 
values of T/\W\ (Centrella, New, Lowe & Brown 2001). 
However, even restricting attention to those configurations 



that are not toriodal and have < 0.27, we find con- 

figurations with masses larger than the maximum mass of 
the corresponding non-rotating star by over 60 %. BSS 
demonstrated that at least some of these models are dy- 
namically stable. Shibata & Uryu (2000) demonstrated 
that binary mergers may result in similarly stable hyper- 
massive stars, when the progenitor masses are not too close 
to the maximum mass. 

To summarize, we find that differential rotation is very 
effective in increasing the maximum allowed mass, espe- 
cially for moderately stiff equations of state. The effect is 
probably large enough to stabilize the remnants of binary 
neutron star merger, which are likely to be differentially 
rotating. Binary neutron star coalescence may therefore 
lead to secularly stable, hypermassive neutron stars. As 
discussed in BSS (see also Shapiro 2000), magnetic braking 
is likely to bring such differentially rotating stars into uni- 
form rotation, which reduces the maximum allowed mass 
and induces a delayed collapse to a Kerr black hole. 
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at the University of Illinois at Urbana-Champaign and 
NSF Grant PHY-0139907 at Bowdoin College. NDL grate- 
fully acknowledges support through an Undergraduate Re- 
search Assistantship from the Department of Physics and 
Astronomy at Bowdoin College. 



APPENDIX 

NUMERICAL RESULTS FOR MAXIMUM MASSES 

We list below in Tables A2 to A9 values for the maximum rest mass increase for uniformly and diff'erentially rotating 

polytropcs. For each polytropic index n we tabulate the differential rotation parameter A~^, the ratio of the central and 
equatorial angular velocity flc/fle (which reduces to (20) in the Newtonian limit, i.e. for soft equations of state), the 
numerically determined fractional rest mass increase {6M/M)num, the ratio of the (relativistic) rotational kinetic energy 
and the gravitational binding energy T/|I^|, the ratio between polar and equatorial radius Rp/Rg, the maximum density 
Pmax, and the estimate {SM/M)est according to equation (17). In these estimates we used the numerically determined 
ratios flc/fle and the central condensations according to (18). For n < 1.25 some of the maximum mass configurations 
are toriodal, i.e. assume the maximum density on a toroid about the center. For these polytropic indices we also include 
the ratio Pc/Pmax- 

All models are computed with the code of CSTl and CST2, using 64 zones both in the radial and angular direction, 
and truncating the Legendre polynomial expansion at £ = 16 (see CSTl for details of the numerical implementation). The 

accuracy of individual stellar models can be tested, for example, by computing a relativistic Virial theorem (Gourgoulhon 
& Bonazzola 1994; Cook, Shapiro & Teukolsky 1996; see also Nozawa et al. 1998 for a comparison of several different 
computational methods). In our analysis, however, the error in the maximum mass and related quantities is dominated 
by the finite step size in the sequences over Rp/Rf, and Pmax, which result in errors typically in the order of a few percent. 
For soft equations of state the mass as a function of central density is a very slowly varying function, making it quite 
difficult to determine the central density of the maximum mass configuration very accurately. The error in the Rp/R^ is 
determined by our stepsize of 0.025. We finally note that highly toroidal configurations depend very sensitively on the 
input parameters, so that those mass increases should only be taken as estimates. 



Table A 2 

n = 0.5: Mo™ = 0.151; p/Pc|tov = 0.375. 



i-l 




{SM/M)num 


T/\W\ 


Rp/Re 


Pmax 


Pmax/Pc 


(5M/M)e.t 


0.0 


1.00 


0.22 


0.15 


0.55 


1.02 


1 


0.38 


0.3 


1.51 


0.41 


0.23 


0.425 


0.73 


1 


0.86 


0.5 


1.93 


0.62 


0.31 


0.2 


0.30 


0.77 


1.40 


0.7 


2.46 


0.46 


0.30 


0.025 


0.29 


0.085 


2.27 


1.0 


3.54 


0.20 


0.27 


0.1 


0.31 


0.29 


4.70 
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Table A3 

n = 0.75: M5°^ = 0.159; p/Pc|tov = 0.28. 





c / 


(5M/M)nnTn 


T/\W\ 


Rp/ Re 






(^M/M)est 


0.0 


1.00 


0.18 


0.11 


0.575 


0.67 


1 


0.28 


0.3 


1.35 


0.27 


0.15 


0.5 


0.64 


1 


0.52 


0.5 


1.92 


0.51 


0.22 


0.4 


0.42 


1 


1.04 


0.7 


2.49 


1.07 


0.30 


0.025 


0.16 


0.021 


1.75 


1.0 


3.48 


0.68 


0.27 


0.025 


0.16 


0.019 


3.42 


1.5 


6.33 


0.19 


0.16 


0.425 


0.27 


0.74 


11.3 








Table A4 












n = 1: Mo™ = 0.180: p/Pc|tov = 0.209. 








{6M/M)num 


T/\W\ 


Rp/ Re 


Pmax 


Pmax/ Pc 


((5M/M)est 


0.0 


1.00 


0.15 


0.083 


0.575 


0.35 


1 


0.21 


0.3 


1.24 


0.20 


0.10 


0.55 


0.33 


1 


0.32 


0.5 


1.65 


0.31 


0.14 


0.475 


0.32 


1 


0.57 


0.7 


2.33 


0.61 


0.21 


0.375 


0.23 


1 


1.14 


0.8 


2.66 


1.12 


0.28 


0.25 


0.083 


0.65 


1.48 


0.85 


2.78 


1.40 


0.29 


0.025 


0.068 


5.0x10^3 


1.60 


1.0 


3.39 


1.22 


0.28 


0.025 


0.075 


4.7x10-3 


2.28 


1.5 


6.33 


0.31 


0.15 


0.475 


0.23 


0.81 


8.37 








Table A5 












n = 1.25: 


M?"^ = 0.216; p/pc\ 


Tov = 0.15. 




i-i 




(5M/M)nu„i 


T/\W\ 


Rp/Re 


Pmax 


Pmax/ pc 


(5M/M)est 


0.0 


1.00 


0.12 


0.063 


0.6 


0.15 


1 


0.15 


0.3 


1.19 


0.16 


0.075 


0.575 


0.15 


1 


0.21 


0.5 


1.51 


0.21 


0.10 


0.525 


0.15 


1 


0.33 


0.7 


1.98 


0.32 


0.13 


0.475 


0.14 


1 


0.58 


1.0 


3.39 


1.78 


0.28 


0.025 


0.037 


1.1x10-3 


1.69 


1.5 


5.26 


0.27 


0.12 


0.575 


0.101 


0.99 


4.10 



Table A6 

n = 1.5: M^"''" = 0.276; p/Pc\toy = 0.115. 



i-l 




{dM/M)nurr 


L T/\W\ 


Rp/Re 


Pmax 


{6M/M),,t 


0.0 


1.00 


0.10 


0.047 


0.625 


0.061 


0.12 


0.3 


1.15 


0.12 


0.055 


0.6 


0.060 


0.15 


0.5 


1.42 


0.16 


0.068 


0.575 


0.059 


0.23 


0.7 


1.81 


0.22 


0.089 


0.525 


0.056 


0.38 


1.0 


2.65 


0.40 


0.137 


0.45 


0.047 


0.81 


1.5 


4.75 


0.25 


0.098 


0.625 


0.050 


2.60 








Table A 7 










n = 2.0: Mo'""" 


— 0.523; p//9c|tov - 


-- 0.063. 










T/\W\ 


Rp/Re 


Pmax X 10^ 


{SM/MU 


0.0 


1.00 


0.067 


0.027 


0.65 


5.1 


0.063 


0.3 


1.12 


0.076 


0.031 


0.625 


5.1 


0.079 


0.5 


1.33 


0.092 


0.036 


0.625 


5.1 


0.11 


0.7 


1.63 


0.12 


0.045 


0.6 


4.9 


0.17 


1.0 


2.28 


0.18 


0.064 


0.55 


4.4 


0.33 


1.5 


4.03 


0.15 


0.059 


0.7 


7.0 


1.03 
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Table A8 

n = 2.5: M^"'^" = 1.25; p/Pc\tov = 0.034. 



i-l 




((5M/M)num 


T/\W\ 


Rp/ Re 


Pmax X 10" 


(5M/M)est 


0.0 


1.00 


0.043 


0.016 


0.675 


1.15 


0.034 


0.3 


1.10 


0.048 


0.017 


0.65 


1.15 


0.041 


0.5 


1.28 


0.056 


0.020 


0.65 


1.1 


0.056 


0.7 


1.54 


0.069 


0.024 


0.625 


1.1 


0.081 


1.0 


2.10 


0.098 


0.034 


0.6 


1.1 


0.15 


1.5 


3.50 


0.102 


0.035 


0.75 


1.05 


0.42 


2.0 


5.44 


0.053 


0.019 


0.875 


1.1 


1.01 



Table A9 

n = 2.9: MJ"*^ = 3.23; p/Pc|tov = 0.021. 



i-l 




(<5M/M)num 


T/\W\ 


Rp/ i?e 


Pmax X 10^ 


(5M/M)est 


0.0 


1.00 


0.028 


0.010 


0.675 


1.4 


0.021 


0.3 


1.09 


0.031 


0.011 


0.675 


1.3 


0.025 


0.5 


1.25 


0.037 


0.013 


0.675 


1.3 


0.033 


0.7 


1.50 


0.044 


0.015 


0.65 


1.3 


0.047 


1.0 


2.02 


0.062 


0.020 


0.625 


1.3 


0.085 


1.5 


3.29 


0.069 


0.022 


0.775 


1.3 


0.227 


2.0 


5.07 


0.035 


0.012 


0.9 
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